Gulf of Maine Sea Surface Temperature Outlook

Gulf of Maine Domain

The spatial extent for the “Gulf of Maine” has been adjusted to meet analysis needs and to match the data properties of different projects. For the CPR analyses the Gulf of Maine region was determined as a bounding box with the following lat/lon bounds lon = c(-70.0, -66.6) & lat = c(42.2, 43.8). For any analyses that rely on groundfish survey data the area representing the Gulf of Maine is a union of individual survey strata. The strata that represent the Gulf of Maine region are strata 24 - 40.

Both of these polygons are meaningful to the data sources they are used with, but this highlights a common issue for maintaining consistent records of temperature and anomalies. What is the source polygon used for any sort of mask, and at what scale are climatologies and anomalies calculated.

####  Data  ####

# Load polygons for mapping
new_england <- ne_states("united states of america") %>% st_as_sf(crs = 4326) 
canada <- ne_states("canada") %>% st_as_sf(crs = 4326)
world <- ne_countries() %>% st_as_sf(crs = 4326)

# File Paths

if(params$env == "local"){
  mills_path <- shared.path(os.use = "unix", group = "Mills Lab", folder = NULL)
  res_path <- shared.path(os.use = "unix", group = "RES Data", folder = NULL)
} else {
  mills_path <-here("Mills Lab/")
  res_path <- here("Res Data/")
}

# Stratum Shapefiles
survey_strata <- read_sf(str_c(res_path, "Shapefiles/BottomTrawlStrata/BTS_Strata.shp"))  %>% 
  clean_names() %>% 
  filter(strata >= 01010 ,
         strata <= 01760,
         strata != 1310,
         strata != 1320,
         strata != 1330,
         strata != 1350,
         strata != 1410,
         strata != 1420,
         strata != 1490) 


# Key to which strata = which regions
strata_key <- list(
  "Georges Bank"          = as.character(13:23),
  "Gulf of Maine"         = as.character(24:40),
  "Southern New England"  = str_pad(as.character(1:12), width = 2, pad = "0", side = "left"),
  "Mid-Atlantic Bight"    = as.character(61:76))


# Assign Areas
strata_assigned <- survey_strata %>% 
  mutate(
    strata = str_pad(strata, width = 5, pad = "0", side = "left"),
    area = case_when(
    str_sub(strata, 3, 4) %in% strata_key$`Georges Bank`         ~ "Georges Bank",
    str_sub(strata, 3, 4) %in% strata_key$`Gulf of Maine`        ~ "Gulf of Maine",
    str_sub(strata, 3, 4) %in% strata_key$`Southern New England` ~ "Southern New England",
    str_sub(strata, 3, 4) %in% strata_key$`Mid-Atlantic Bight`   ~ "Mid-Atlantic Bight",
    TRUE ~ "Outside Major Study Areas"
  )) 



# Focus on Gulf of Maine Zones
gm  <- strata_assigned %>% filter(area == "Gulf of Maine")

# Re-project
utm_19      <- "+proj=utm +zone=19 +ellps=WGS84 +datum=WGS84 +units=m +no_defs "
gm_utm      <- gm %>% st_transform(crs = utm_19)
canada_utm  <- canada %>% st_transform(crs = utm_19)
usa_utm     <- new_england %>% st_transform(crs = utm_19)

# Crop bounds
crop_x <- st_bbox(gm_utm)[c(1,3)] 
crop_y <- st_bbox(gm_utm)[c(2,4)]

# Adjust bbox (units are meters)
crop_x <- crop_x + c(-30000, 20000)
crop_y <- crop_y + c(-20000, 20000)


# Map the survey stratum
survey_p <- ggplot() +
  geom_sf(data = usa_utm) +
  geom_sf(data = canada_utm) +
  geom_sf(data = gm_utm, aes(fill = area), alpha = 0.4, show.legend = F) +
  coord_sf(xlim = crop_x, ylim = crop_y, expand = T) +
  labs(title = "Groundfish Survey", subtitle = "Gulf of Maine Region")



####  Cpr Survey Region

# From: adamkemberling/continuous_plankton_recorder 
# CPR Extent
cpr_extent <- tribble(
  ~"lon",      ~"lat",
  -70.000000,   42.200000, 
  -66.600000,   42.200000, 
  -66.600000,   43.800000, 
  -70.000000,   43.800000, 
  -70.000000,   42.200000
) %>% 
  as.matrix() %>% 
  list() %>% 
  st_polygon()

# as simple feature collection
gm_cpr <- st_sf(area = "CPR Gulf of Maine", st_sfc(cpr_extent), crs = 4326)
cpr_utm <- st_transform(gm_cpr, crs = utm_19)


# Plot with the rest
cpr_p <- ggplot() +
  geom_sf(data = usa_utm) +
  geom_sf(data = canada_utm) +
  geom_sf(data = cpr_utm, aes(fill = area), alpha = 0.4, show.legend = F) +
  coord_sf(xlim = crop_x, ylim = crop_y, expand = T) +
  labs(title = "CPR Survey", subtitle = "Gulf of Maine Region")





####  Ecological Production units  ####
epu_sf <- read_sf(str_c(res_path, "Shapefiles/EPU/EPU_extended.shp"), crs = 4326)
epu_utm <- st_transform(epu_sf, crs = utm_19)
epu_utm <- filter(epu_utm, EPU == "GOM")

epu_p <- ggplot() +
  geom_sf(data = usa_utm) +
  geom_sf(data = canada_utm) +
  geom_sf(data = epu_utm, aes(fill = EPU), alpha = 0.4, show.legend = F) +
  coord_sf(xlim = crop_x, ylim = crop_y, expand = T) +
  labs(title = "Ecological Production Units", subtitle = "Gulf of Maine Region")



####  Physio Regions  ####
nelme_sf <- read_sf(str_c(res_path, "Shapefiles/NELME_regions/GoM_sf.shp"), crs = 4326)
nelme_utm <- st_transform(nelme_sf, crs = utm_19)

nelme_p <- ggplot() +
  geom_sf(data = usa_utm) +
  geom_sf(data = canada_utm) +
  geom_sf(data = nelme_utm, aes(fill = "GOM"), alpha = 0.4, show.legend = F) +
  coord_sf(xlim = crop_x, ylim = crop_y, expand = T) +
  labs(title = "NELME", subtitle = "Gulf of Maine Region")

# patchwork
(survey_p | cpr_p) / (epu_p | nelme_p) & theme(axis.text = element_text(size = 8))

Accessing OISST on Box

For many of our most-used areas climatologies and anomalies have been processed using a handful of jupyter notebook workflows to take advantage of {xarray}. The following sections will highlight what those workflows are, what they detail, and how they can be updated using the Gulf of Maine as an example.

The root folder for the following OISST products is ~/Box/NSF OKN Demo Data/oisst.

Masked Timeseries

A masked timeseries is the most basic building block for any of these areas. For any desired area (represented by a spatial polygon) a timeseries table of observed sea surface temperature is returned that contains the mean sea surface temperature over that area at each time step. Rather than repeat this process again for the climatology table and the anomalies these three values and a few additional fields are stored in the same tables.

Pre-processed tables have been placed on box as part of the NSF C-ACCEL project under the folder: ~/Box/NSF OKN Demo Data/oisst/likelihood_timeseries

The naming convention in the folder is: likelihood_timeseries/nmfs_trawl_regions/OISST_v2_anom_ + region_name + _likelihood_ts.csv

if (params$env == "local"){
gom_path <- "~/Box/NSF OKN Demo Data/oisst/likelihood_timeseries/nmfs_trawl_regions/OISSTv2_anom_gulf_of_maine_likelihood_ts.csv"
} else{
  gom_path <- "NSF OKN Demo Data/oisst/likelihood_timeseries/nmfs_trawl_regions/OISSTv2_anom_gulf_of_maine_likelihood_ts.csv"
}
  
gom_ts <- read_csv(gom_path)

head(gom_ts) %>% kable() %>% kable_styling()
time modified_ordinal_day sst clim_sd sst_clim sst_anom log_lik
1981-09-01 244 15.56767 1.712336 16.44248 -0.8748112 1.587300
1981-09-02 245 15.63219 1.753561 16.29203 -0.6598396 1.551383
1981-09-03 246 15.43664 1.742908 16.24344 -0.8067951 1.581632
1981-09-04 247 14.85089 1.714197 16.18611 -1.3352232 1.761242
1981-09-05 248 14.67466 1.683289 16.14633 -1.4716759 1.821876
1981-09-06 249 14.85329 1.613693 16.00317 -1.1498871 1.651350

Climatologies are currently set up to calculate daily averages on a modified julian day, such that every March 1st and all days after fall on the same day, regardless of whether it is a leap year or not. In these tables sst is the mean temperature observed for that date averaged across all cells within the area. sst_clim & clim_sd are the climate means and standard deviations for a 1982-2011 climatology. sst_anom is the daily observed minus the climate mean, and log_lik is the log likelihood of getting that value given a normal distribution with mean of sst_clim and sd of clim_sd.

Marine Heatwaves

The {heatwaveR} package provides a relatively quick way of working with tabular data to calculate a seasonal climate mean, as well as heatwaves and coldwaves at a desired threshold. These functions can be used to get the climatology using the standard day of year if desired, and track the number and length of heatwave events. Heatwave events follow the definition of Hobday et al. 2016

The heatwave history for the Gulf of Maine trawl region (~from the NELME regions shapefile~) is as follows:

# Wrapper function to do heatwaves and coldwaves simultaneously at 90%
pull_events <- function(temperature_timeseries) {
    
    # Pull the two column dataframe for mhw estimation
    test_ts <- data.frame(t = temperature_timeseries$time, 
                          temp = temperature_timeseries$sst)
    
    # Detect the events in a time series
    ts  <- ts2clm(data = test_ts, climatologyPeriod = c("1982-01-01", "2011-12-31"))
    
    #heatwaves
    mhw <- detect_event(ts)                         
    
    # prep heatwave data
    mhw_out <- mhw$climatology %>% 
        mutate(sst_anom = temp - seas) %>% 
        select(time = t,
               sst = temp,
               seas,
               sst_anom,
               mhw_thresh = thresh,
               mhw_event = event,
               mhw_event_no = event_no)

    # 2. Detect cold spells
    ts <- ts2clm(data = test_ts, 
                 climatologyPeriod = c("1982-01-01", "2011-12-31"), 
                 pctile = 10)
    mcs <- detect_event(ts, coldSpells = TRUE) #coldSpells = TRUE flips boolean to < thresh
    
    # prep cold spell data
    mcs_out <- mcs$climatology %>%
        select(time = t,
               mcs_thresh = thresh,
               mcs_event = event,
               mcs_event_no = event_no)
    

    # join heatwaves to coldwaves
    hot_and_cold <- left_join(mhw_out, mcs_out, by = "time")
    
    
    
    # 3. Data formatting for plotting, 
    # adds columns to plot hw and cs seperately
    events_out <- hot_and_cold %>% 
        mutate(status = ifelse(mhw_event == TRUE, "Marine Heatwave Event", "Sea Surface Temperature"),
               status = ifelse(mcs_event == TRUE, "Marine Cold Spell Event", status),
               hwe = ifelse(mhw_event == TRUE, sst, NA),
               cse = ifelse(mcs_event == TRUE, sst, NA),
               nonevent = ifelse(mhw_event == FALSE & mcs_event == FALSE, sst, NA)) 
    
    # Close the gaps between a mhw event and sst (might not need if full line for temp exists)
    events_out <- events_out %>% 
        mutate(hwe = ifelse(is.na(hwe) & is.na(lag(hwe)) == FALSE, sst, hwe),
               cse = ifelse(is.na(cse) & is.na(lag(cse)) == FALSE, sst, cse))
    
    
    return(events_out)
}



# Use function
#
gom_heatwaves <- pull_events(gom_ts)

Heatwave timelines can be interactively plotted using {ggplot2} for static plots or {plotly} for interactive content.

Interactive Plot

# Helper function
plotly_mhw_plots <- function(data){
    
    # How to get rgb() from colors
    plot_cols <- list(
        "gray20"    = col2rgb(col = "gray20"),
        "gray40"    = col2rgb(col = "gray40"),
        "royalblue" = col2rgb(col = "royalblue"),
        "darkred"   = col2rgb(col = "darkred"),
        "royalblue" = col2rgb(col = "royalblue"))
    
    
    # Building the plot
    fig <- plot_ly(data, x = ~time) 
    
    # Sea Surface Temperature
    fig <- fig %>% add_trace(y = ~sst, 
                             name = 'Sea Surface Temperature',
                             mode = 'lines', 
                             type = "scatter",
                             line = list(color = "rgb(65, 105, 225)", 
                                         width = 2)) 
    # Heatwave Threshold
    fig <- fig %>% add_trace(y = ~mhw_thresh, 
                             name = 'MHW Threshold ', 
                             mode = 'lines', 
                             type = "scatter",
                             line = list(color = "rgb(51, 51, 51)", 
                                         width = 1, 
                                         dash = 'dot')) 
    # Seasonal Climatology
    fig <- fig %>% add_trace(y = ~seas, 
                             name = 'Seasonal Climatology ', 
                             mode = 'lines', 
                             type = "scatter",
                             line = list(color = "rgb(102, 102, 102)", 
                                         width = 3, 
                                         dash = 'dash')) 
    # Marine Cold Spell Threshold
    fig <- fig %>% add_trace(y = ~mcs_thresh, 
                             name = 'MCS Threshold ', 
                             mode = 'lines', 
                             type = "scatter",
                             line = list(color = "rgb(51, 51, 51)", 
                                         width = 1, 
                                         dash = 'dot')) 
    # Heatwave Event
    fig <- fig %>% add_trace(y = ~hwe, 
                             name = 'Marine Heatwave Event ', 
                             mode = 'lines', 
                             type = "scatter",
                             line = list(color = "rgb(139, 0, 0)", 
                                         width = 2)) 
   
    # Cold Spell Event
    fig <- fig %>% add_trace(y = ~cse, 
                             name = 'Marine Cold Spell Event ', 
                             mode = 'lines', 
                             type = "scatter",
                             line = list(color = "rgb(65, 105, 225)", 
                                         width = 2)) 
    
    
    fig <- fig %>% layout(xaxis = list(title = ""),
                          yaxis = list (title = "Temperature (degrees C)"))
    return(fig)
}



# Plot the last year
last_year <- Sys.Date() - 365
gom_heatwaves %>% 
  filter(time >= last_year) %>% 
  plotly_mhw_plots()

Static Plot

color_vals <- c(
  "Sea Surface Temperature" = "royalblue",
  "Heatwave Event"          = "darkred",
  "Cold Spell Event"        = "lightblue",
  "MHW Threshold"           = "gray30",
  "MCS Threshold"           = "gray30",
  "Seasonal Climatology"    = "gray30")


ylab <- expression("Sea Surface Temperature"~degree~C)

gom_heatwaves %>% 
  filter(time >= last_year) %>% 
  ggplot(aes(x = time)) +
    geom_segment(aes(x = time, xend = time, y = seas, yend = sst), color = "royalblue", alpha = 0.25) +
    geom_segment(aes(x = time, xend = time, y = mhw_thresh, yend = hwe), color = "darkred", alpha = 0.25) +
    geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
    geom_line(aes(y = hwe, color = "Heatwave Event")) +
    geom_line(aes(y = cse, color = "Cold Spell Event")) +
    geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 2, size = .5) +
    geom_line(aes(y = mcs_thresh, color = "MCS Threshold"), lty = 2, size = .5) +
    geom_line(aes(y = seas, color = "Seasonal Climatology"), size = .5) +
    scale_color_manual(values = color_vals) +
    theme(legend.title = element_blank()) +
    labs(x = "", y = ylab)

Heatwave Events

# prep the legend title
guide_lab <- expression("Sea Surface Temperature Anomaly"~degree~C)

# get new axis dimensions, y = year, x = day within year
# use flate_date so that they don't stair step
base_date <- as.Date("2000-01-01")
grid_data <- gom_heatwaves %>% 
  mutate(year = year(time),
         yday = yday(time),
         flat_date = as.Date(yday-1, origin = base_date))


# assemble plot
ggplot(grid_data, aes(x = flat_date, y = year)) +
  # background box fill for missing dates
  geom_rect(xmin = base_date, xmax = base_date + 365, 
            ymin = min(grid_data$year) - .5, ymax = max(grid_data$year) + .5, 
            fill = "gray75", color = "transparent") +
  # tile for sst colors
  geom_tile(aes(fill = sst_anom)) +
  # points for heatwave events
  geom_point(data = filter(grid_data, mhw_event == TRUE),
             aes(x = flat_date, y = year), size = .25)  +
  scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
  scale_y_continuous(limits = c(1980.5, 2020.5), expand = c(0,0)) +
  labs(x = "", y = "") +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
  #5 inches is default rmarkdown height
  guides("fill" = guide_colorbar(title = guide_lab, 
                                 title.position = "right", 
                                 title.hjust = 0.5,
                                 barheight = unit(4.8, "inches"), 
                                 frame.colour = "black", 
                                 ticks.colour = "black")) +  
  theme_classic() +
  theme(legend.title = element_text(angle = 90))

Anomaly Maps

For demonstration purposes I will use 2020 anomaly data.

# Access information to netcdf on box
okn_path <- shared.path(group = "NSF OKN", folder = "oisst/annual_anomalies/")
nc_year <- "2020"
anom_path <- str_c(okn_path, "daily_anoms_", nc_year, ".nc")

# Load 2020 as stack
anoms_2020 <- stack(anom_path)

Mean Global

The following code will subset the anomalies for July and plot them:

# Get the mean temperature anomalies for July
july_dates <- which(str_sub(names(anoms_2020), 7, 8) == "07")
july_avg <- mean(anoms_2020[[july_dates]])
#plot(july_avg, main = "Mean Temperature Anomaly - July 2020")


# Convert wgs84 to stars
july_st <- st_as_stars(rotate(july_avg))



####  NOTE: having trouble projecting to either robinson or equal earth projections for some resson  ####
# Solved : Tranformation vs. Warping https://r-spatial.github.io/stars/articles/stars5.html

# robinson projection
robinson_proj <- "+proj=robin"

# equal earth projection - not working, don't think sf has functionality for it yet
eqearth_proj <- "+proj=eqearth"

# warp to equal robinson or earth projection
robinson_grid <- july_st %>% 
  st_transform(robinson_proj) %>% 
  st_bbox() %>%
  st_as_stars()
july_robinson <- july_st %>% 
  st_warp(robinson_grid) 



# Plot
ggplot() +
  geom_stars(data = july_st) +
  geom_sf(data = world, fill = "gray90") +
  scale_fill_gradient2(low = "blue",
                       mid = "white",
                       high = "red") +
  map_theme +
  coord_sf(expand = FALSE) +
  guides("fill" = guide_colorbar(title = "Average Sea Surface Temperature Anomaly",
                                 title.position = "top", 
                                 title.hjust = 0.5,
                                 barwidth = unit(4, "in"), 
                                 frame.colour = "black", 
                                 ticks.colour = "black"))

# ####  NOTE: having trouble projecting to either robinson or equal earth projections for some resson  ####
# # Solved : Tranformation vs. Warping https://r-spatial.github.io/stars/articles/stars5.html
# 
# # robinson projection
# robinson_proj <- "+proj=robin"
# 
# # equal earth projection - not working, don't think sf has functionality for it yet
# eqearth_proj <- "+proj=eqearth"
# 
# # warp to equal robinson or earth projection
# robinson_grid <- july_st %>% 
#   st_transform(robinson_proj) %>% 
#   st_bbox() %>%
#   st_as_stars()
# july_robinson <- july_st %>% 
#   st_warp(robinson_grid) 
# 
# 
# 
# # Plot - issues with land masses not matching raster
# ggplot() +
#   geom_stars(data = july_robinson) +
#   geom_sf(data = st_transform(world, robinson_proj), fill = "gray90") +
#   scale_fill_gradient2(low = "blue",
#                        mid = "white",
#                        high = "red") +
#   map_theme +
#   coord_sf(expand = FALSE) +
#   guides("fill" = guide_colorbar(title = "Average Sea Surface Temperature Anomaly",
#                                  title.position = "top", 
#                                  title.hjust = 0.5,
#                                  barwidth = unit(4, "in"), 
#                                  frame.colour = "black", 
#                                  ticks.colour = "black"))

GoM

Using the July average sst anomaly we can then clip down to the Gulf of Maine for an aerial view of the region:

# Clip Raster - Convert to stars
gom_ras <- crop(rotate(july_avg), extent(-77, -60, 35, 46))
gom_st <- st_as_stars(gom_ras)

# Crop bounds 
crop_x <- st_bbox(gom_st)[c(1,3)] 
crop_y <- st_bbox(gom_st)[c(2,4)]

# Plot Gulf of Maine - wgs84
ggplot() +
  geom_stars(data = gom_st) +
  geom_sf(data = new_england, fill = "gray90") +
  geom_sf(data = canada, fill = "gray90") +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", na.value = "transparent") +
  map_theme +
  coord_sf(xlim = crop_x, ylim = crop_y, expand = F) +
  guides("fill" = guide_colorbar(
    title = "Average Sea Surface Temperature Anomaly",
    title.position = "top",
    title.hjust = 0.5,
    barwidth = unit(4, "in"),
    frame.colour = "black",
    ticks.colour = "black"))

####  Albers Projection - Working  ####
# albers equal area: centered on -70 degrees
# The settings for lat_1 and lat_2 are the locations at which the cone intersects the earth, so distortion is minimized at those latitudes
alb_70 <- "+proj=aea +lat_1=30 +lat_2=50 +lon_0=-70"


# # Projecting a raster - results in curvilinear grid - not plottable in ggplot
# gom_alb <- st_transform(gom_st, st_crs(alb_70))
# plot(gom_alb)


# warping instead of transforming:
# create regular grid in new coordinate reference system
# then warp values over
newgrid <- gom_st %>% 
  st_transform(alb_70) %>% 
  st_bbox() %>%
  st_as_stars()
gom_alb <- gom_st %>% 
  st_warp(newgrid) 

# transform polygons because we will use them again
ne_alb <- st_transform(new_england, alb_70)
can_alb <- st_transform(canada, alb_70)
world_alb <- st_transform(world, alb_70)

# Crop bounds using albers projection
crop_x <- st_bbox(gom_alb)[c(1,3)] 
crop_y <- st_bbox(gom_alb)[c(2,4)]


gom_albers_plot <- ggplot() +
  geom_stars(data = gom_alb) +
  geom_sf(data = ne_alb, fill = "gray90") +
  geom_sf(data = can_alb, fill = "gray90") +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", na.value = "transparent") +
  map_theme +
  coord_sf(crs = alb_70, xlim = crop_x, ylim = crop_y, expand = T) +
  guides("fill" = guide_colorbar(
    title = "Average Sea Surface Temperature Anomaly",
    title.position = "top", 
    title.hjust = 0.5,
    barwidth = unit(4, "in"), 
    frame.colour = "black", 
    ticks.colour = "black"))

NW Atlantic

I was not able to locate a shapefile for the Northwest Atlantic Region so I’m going to eyeball it for now off the image in the word doc. Good practice for making masks or shapes programmatically.

Andy’s changes:
> I’d suggest we move the southern extent down to at least Long Island, and we may want to move the eastern bound a bit further to pick up the southern part of E Greenland (~over to 40 W long).

nw_extent <- tribble(
  ~"lon",  ~"lat",
  -73,  40.5, 
  -40,  40.5, 
  -40,  70, 
  -73,  70, 
  -73,  40.5
) %>% as.matrix() %>% 
  list() %>% 
  st_polygon()

# create simple feature collection
nw_region <- st_sf(area = "Northwest Atlantic Testing", st_sfc(nw_extent), crs = 4326)
greenland <- ne_states(country = "greenland") %>% st_as_sfc(crs = 4326)

# new sterographic projection for this area, centered using lon_0
stereographic_north <- "+proj=stere +lat_0=90 +lat_ts=75 +lon_0=-57"
nw_north            <- st_transform(nw_region, crs = stereographic_north)
canada_north        <- st_transform(canada, stereographic_north)
newengland_north    <- st_transform(new_england, stereographic_north)
greenland_north     <- st_transform(greenland, stereographic_north)

# coord_sf Crop bounds in projection units
crop_x <- st_bbox(nw_north)[c(1,3)] 
crop_y <- st_bbox(nw_north)[c(2,4)]



ggplot() +
  geom_sf(data = nw_north, aes(fill = area), alpha = 0.3) +
  geom_sf(data = newengland_north, fill = "gray90") +
  geom_sf(data = canada_north, fill = "gray90") +
  geom_sf(data = greenland_north, fill = "gray90") +
  coord_sf(crs = stereographic_north, xlim = crop_x, ylim = crop_y, expand = T) +
  guides("fill" = guide_legend(title = "")) +
  theme(
    panel.border = element_rect(color = "black", fill = NA),
    plot.background = element_rect(color = "transparent", fill = "transparent"),
    axis.title.x = element_blank(), # turn off titles
    axis.title.y = element_blank(),
    legend.position = "bottom", 
    legend.title.align = 0.5)

NW Atlantic OISST

 

A work by Adam A. Kemberling

Akemberling@gmri.org